Multiple Numerical Assembly

Sections

Abstract

In this document, I want to provide a write-up of how this model is similar to and different from the versions previously implemented. This model is meant to provide assembly of multiple sites at the same time, that may, or may not, be connected in some fashion. Along the way, I will mention some of the input and output format that I am expect as documentation. At the end, I will be presenting some preliminary results. I especially want to use this as a vehicle to think about how to analyse those results.

Functions

library(dplyr)     # Data Manipulation

Warning: package ‘dplyr’ was built under R version 4.1.1

Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

library(ggplot2)   # 2-D Plot

Warning: package ‘ggplot2’ was built under R version 4.1.1

library(plotly)    # 3-D Plot

Warning: package ‘plotly’ was built under R version 4.1.1
Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

Attaching package: ‘plotly’

The following object is masked from ‘package:ggplot2’:

    last_plot

The following object is masked from ‘package:stats’:

    filter

The following object is masked from ‘package:graphics’:

    layout

library(ggfortify) # used for biplots of PCAs

Warning: package ‘ggfortify’ was built under R version 4.1.1

library(vegan)     # Ecological analysis mega-package

Warning: package ‘vegan’ was built under R version 4.1.1
Loading required package: permute
Warning: package ‘permute’ was built under R version 4.1.1
Loading required package: lattice
This is vegan 2.5-7

library(RMTRCode2)

# https://stackoverflow.com/a/7172832
ifrm <- function(obj, env = globalenv()) {
    obj <- deparse(substitute(obj))
    if(exists(obj, envir = env)) {
        rm(list = obj, envir = env)
    }
}

Results, No Spatial Structure

First, we load up the preliminary results. In this case, we are loading a system in which 34 basal species and 66 consumer species form the pool for 10 unconnected environments. The pool and interaction matrix were assembled with the default parameters from Law and Morton’s 1996 work.

load(file.path(
  "..", "experiments", "MNA-FirstAttempt-Result-Env10-None.RData")
)
toRemove <- result$Abundance <= result$Parameters$EliminationThreshold
result$Abundance[toRemove] <- 0
ifrm(toRemove)

Events

In total, 9320 events were used in these environments, with the species and environment invasion both randomly assigned. The number of arrival and extinction events were controlled to both be half of this number. We chose this number due to the coupon collecting problem. In particular, we use the result that the probability of encountering each species is bounded: \[\text{Pr}(\text{Draws} < n \log_{e} n + c n) \rightarrow \exp(-\exp(-c)) \text{ as } n \rightarrow \infty\] where \(n\) is the number of species and \(c\) is a constant. For our purposes, we choose \(c = 5\) so that we have a probability of about \(99.3\%\) of seeing each species in each environment. In practice, we failed to observe 2 species-environment combinations. Notably, nearly every species had at least one successful invasion; 17 did not.

The initial abundance was set to be 4000 times the elimination threshold, in line with work on minimum viable populations (Traill et al. 2007). The elimination threshold is admittedly more arbitrary, since it sets an effective individual-area relationship. For this calculation, I set it to \(10^{-4}\), in line with our previous calculations. This is large enough to avoid numerical difficulties from precision, while being low enough to represent a decent sized region.

Since each population is assembling simultaneously, I chose to use exponential waiting times for the inter-species arrival and extinction times. Note that these rates are shared between species and environments, but arrival and extinction are fully independent of each other. Species and environment affected in each event was chosen uniformly at random. The question then is how to set the rate.

To set the rate in this case, I chose to set it to the largest eigenvalue magnitude of the per-environment interaction matrices. This magnitude corresponds to the strongest response we can see from the interaction matrix and, if the interaction matrix is a good approximation for the Jacobian around a stable fixed point (which is not guaranteed), indicates the characteristic time scale of the decay to equilibrium. Hence, (overall) arrival rates and (overall) extinction rates should happen on the same timescale as the (largest) dynamics in the system. Since there are 10 environments, we should then expect that 10 characteristic time scales, on average, should occur in between arrival events in the same environment.

Abundance

With 10 environments, it is probably not helpful to check 10 individual abundance curves, but looking at the first one might be helpful.

LawMorton1996_PlotAbundance(result$Abundance[seq(from = 1, 
                                                 to = nrow(result$Abundance), 
                                                 by = 10), c(1, 2:101)]) -> obj;

obj + ggplot2::scale_y_log10() + ggplot2::guides(color = FALSE)

Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
Warning: Transformation introduced infinite values in continuous y-axis

Every vertical line is a species introduction or extinction by neutral dynamics.

Alpha Diversity plots

Intro

Perhaps more intriguing might be some sense of the biodiversity that we have in each system. We break the abundance results up by environment, then calculate the number of non-zero abundance curves at each time point. We also calculate the Shannon entropy (reminder: higher entropy means more uncertainty which means a flatter distribution).

Diversity <- lapply(
  1:result$NumEnvironments,
  function(i, abund, numSpecies) {
    time <- abund[, 1]
    env <- abund[, 1 + 1:numSpecies + numSpecies * (i - 1)]
    richness <- rowSums(env != 0)
    abundSum <- rowSums(env)
    entropy <- env / abundSum
    entropy <- - apply(
      entropy, MARGIN = 1,
      FUN = function(x) {
        sum(ifelse(x != 0, x * log(x), 0))
      })
    species <- apply(
      env, MARGIN = 1,
      FUN = function(x) {
        toString(which(x > 0))
      }
    )
    data.frame(Time = time, 
               Richness = richness, 
               Entropy = entropy,
               Species = species,
               Environment = i)
  },
  abund = result$Abundance,
  numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)


Diversity <- dplyr::bind_rows(Diversity)
Diversity <- Diversity %>% dplyr::mutate(
  Evenness = Entropy / log(Richness)
)

Richness

ggplot2::ggplot(
  Diversity, 
  ggplot2::aes(
    x = Time,
    y = Richness,
    color = factor(Environment),
    alpha = ifelse(Environment == 1, 1, 0.3)
  )
) + ggplot2::geom_line(
) + ggplot2::geom_line(
  data = Diversity %>% dplyr::group_by(
    Time
    ) %>% dplyr::summarise(
      Richness = mean(Richness)
    ),
  mapping = ggplot2::aes(
    x = Time,
    y = Richness
  ),
  color = "black",
  inherit.aes = FALSE
) + ggplot2::guides(
  alpha = "none"
) + ggplot2::scale_color_discrete(
  "Environment"
)

So richness hovers around a similar regime throughout the majority of the simulation. Note in this plot that we have emphasised one environmental curve and superimposed the mean in black. We do manage to reach heights of 11 species in one environment, but these heights are shortlived. Instead, we seem to observe a (time and environment averaged) value:

mean(Diversity$Richness)

[1] 4.458893

(If we consider the first 10,000 time units as burn-in, we instead see a value of

mean(Diversity$Richness[Diversity$Time > 10000])

[1] 4.587217

Entropy

ggplot2::ggplot(
  Diversity, 
  ggplot2::aes(
    x = Time,
    y = Entropy,
    color = factor(Environment),
    alpha = ifelse(Environment == 1, 1, 0.3)
  )
) + ggplot2::geom_line(
) + ggplot2::geom_line(
  data = Diversity %>% dplyr::group_by(
    Time
    ) %>% dplyr::summarise(
      Entropy = mean(Entropy)
    ),
  mapping = ggplot2::aes(
    x = Time,
    y = Entropy
  ),
  color = "black",
  inherit.aes = FALSE
) + ggplot2::guides(
  alpha = "none"
) + ggplot2::scale_color_discrete(
  "Environment"
)

Warning: Removed 8043 row(s) containing missing values (geom_path).
Warning: Removed 2297 row(s) containing missing values (geom_path).

Entropy has a similar, but highly erratic, behaviour. If one tries to follow one of the entropy curves, then one sees that they have fairly substantial periods of almost smooth behaviour followed by suddenly very noisy behaviour, and noise seems to be the dominant mode if one tries to examine the low alpha environment in the background. There are some easy to make predictions. Extinctions reduce the entropy in the system, as you become more certain about what remains. Analogously, arrivals increase the entropy. We can probably better see the relationship beyond these principles by plotting entropy against richness and connecting observations by environment and time.

Entropy-Richness

# ggplot2::ggplot(
#   Diversity %>% dplyr::filter(Environment < 4), 
#   ggplot2::aes(
#     x = Richness,
#     y = Entropy,
#     group = Environment,
#     color = Time
#   )
# ) + ggplot2::geom_path(
# ) + ggplot2::guides(
#   alpha = "none"
# )

plotly::plot_ly(data = Diversity %>% dplyr::filter(Environment < 2),
                x = ~Richness, y = ~Entropy, z = ~Time, type = "scatter3d",
                mode = "lines", opacity = 1, line = list(color = ~Time))

It seems quite hard to tell, but there does not appear to be any particular orientation (clockwise, counter-clockwise) or similar pattern here.

Evenness

ggplot2::ggplot(
  Diversity, 
  ggplot2::aes(
    x = Time,
    y = Evenness,
    color = factor(Environment),
    alpha = ifelse(Environment == 1, 1, 0.3)
  )
) + ggplot2::geom_line(
) + ggplot2::geom_line(
  data = Diversity %>% dplyr::group_by(
    Time
    ) %>% dplyr::summarise(
      Evenness = mean(Evenness)
    ),
  mapping = ggplot2::aes(
    x = Time,
    y = Evenness
  ),
  color = "black",
  inherit.aes = FALSE
) + ggplot2::guides(
  alpha = "none"
) + ggplot2::scale_color_discrete(
  "Environment"
)

Warning: Removed 10986 row(s) containing missing values (geom_path).
Warning: Removed 2831 row(s) containing missing values (geom_path).

Evenness helps highlight that this is a high variance process but with a relatively constrained mean.

Environment Diversity

We can, of course, flip the idea on its head. Instead of examining the diversity of species within environments, we can look at the diversity of environments occupied by species. Since very few species end up occupying environments, we just look at richness. Unfortunately, this is quite a memory exhaustive task.

EnvDiversity <- lapply(
    1:((ncol(result$Abundance) - 1) / result$NumEnvironments),
    function(i, abund, numSpecies) {
        time <- abund[, 1]
        env <- abund[, 1 + i + numSpecies * (1:result$NumEnvironments - 1)]
        richness <- rowSums(env != 0)
        abundSum <- rowSums(env)
        environments <- apply(
            env, MARGIN = 1,
            FUN = function(x) {
                toString(which(x > 0))
            }
        )
        data.frame(Time = time, 
                   Richness = richness, 
                   Abundance = abundSum,
                   Species = i,
                   Environments = environments)
    },
    abund = result$Abundance,
    numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)

EnvDiversity <- dplyr::bind_rows(EnvDiversity)

ggplot2::ggplot(
  EnvDiversity %>% dplyr::filter(Richness > 1), 
  aes(x = Time, y = Richness, color = Species)
) + geom_point(
  alpha = 0.01, size = 3
) + guides(
  color = "none"
)

One immediately interesting trend here is that very few species are present across more than 5 environments at a given time. Indeed, only

unique(EnvDiversity$Species[EnvDiversity$Richness > 5])

[1] 12 14 22 30

are ever present in more than 5 environments at once. We can also examine how long these periods occur for by species by tabulation. We block times so that entries are all of the same unit length, and round the average richness during the given time unit.

with(EnvDiversity %>% mutate(
  Time = floor(Time)
  ) %>% group_by(
    Time, Species
    ) %>% summarise(
      Richness = round(mean(Richness)),
      .groups = "drop"
      ),
     table(Species, Richness))

       Richness
Species     0     1     2     3     4     5     6     7
    1   65970 14384   439     0     0     0     0     0
    2   34776 31378 11352  3287     0     0     0     0
    3   63206 13000  1767  2820     0     0     0     0
    4   46367 21078 11443  1617   288     0     0     0
    5   30648 31553 16542  1341   709     0     0     0
    6   72061  8732     0     0     0     0     0     0
    7   75522  5271     0     0     0     0     0     0
    8   64759 15508   526     0     0     0     0     0
    9   77919  2874     0     0     0     0     0     0
    10  17817 38679 20192  4105     0     0     0     0
    11  52666 24892  3235     0     0     0     0     0
    12   9669  2676  9781 26641 18597  9922  2294  1213
    13  64640 15703   450     0     0     0     0     0
    14  12094 21314 14556 16445 15148   931   305     0
    15  68564 10732  1497     0     0     0     0     0
    16  69690 10903   200     0     0     0     0     0
    17  30636 33714 15505   938     0     0     0     0
    18  69214 11322   257     0     0     0     0     0
    19  62460 17599   734     0     0     0     0     0
    20  74867  3685  1849   392     0     0     0     0
    21  47906 25702  4185  3000     0     0     0     0
    22   3180 13477 20151 21766 16740  5305   174     0
    23  65137 14940   716     0     0     0     0     0
    24  65792 12313  2688     0     0     0     0     0
    25  66085 11109  3599     0     0     0     0     0
    26  41654 19402 19201   536     0     0     0     0
    27  48808 20447 10279  1259     0     0     0     0
    28  68983 11810     0     0     0     0     0     0
    29  34631 26993 17044  2125     0     0     0     0
    30   5561 19065 20682 11364 11597  7963  3301  1260
    31  65019 15665   109     0     0     0     0     0
    32  19327 26015 29600  5851     0     0     0     0
    33  15128 23405 24345 10059  7856     0     0     0
    34  76432  4361     0     0     0     0     0     0
    35  73829  6080   884     0     0     0     0     0
    36  71050  9743     0     0     0     0     0     0
    37  38199 22369 15663  4486    76     0     0     0
    38  28752 35327 16714     0     0     0     0     0
    39  77008  3785     0     0     0     0     0     0
    40  49953 27129  3711     0     0     0     0     0
    41  79838   955     0     0     0     0     0     0
    42  78312  2481     0     0     0     0     0     0
    43  32813 23892 17056  6661   371     0     0     0
    44  39700 23980 15554  1559     0     0     0     0
    45  80793     0     0     0     0     0     0     0
    46  48461 27669  4663     0     0     0     0     0
    47  12053 36783 14826 15182  1949     0     0     0
    48  63151 16980   662     0     0     0     0     0
    49  73125  7668     0     0     0     0     0     0
    50  79860   933     0     0     0     0     0     0
    51  80255   538     0     0     0     0     0     0
    52  73946  6847     0     0     0     0     0     0
    53  77581  3212     0     0     0     0     0     0
    54  76399  4394     0     0     0     0     0     0
    55  55217 23951  1625     0     0     0     0     0
    56  72564  8229     0     0     0     0     0     0
    57  57443 22644   706     0     0     0     0     0
    58  78988  1805     0     0     0     0     0     0
    59  35315 37942  7536     0     0     0     0     0
    60  67881 12912     0     0     0     0     0     0
    61  52226 24749  3818     0     0     0     0     0
    62  62517 16663  1613     0     0     0     0     0
    63  46937 31441  2415     0     0     0     0     0
    64  68974 11213   606     0     0     0     0     0
    65  80793     0     0     0     0     0     0     0
    66  67130 12561  1102     0     0     0     0     0
    67  37887 27036 11834  4036     0     0     0     0
    68  51139 27047  2205   402     0     0     0     0
    69  61791 16492  2510     0     0     0     0     0
    70  73231  7354   208     0     0     0     0     0
    71  22117 22512 20026 14172  1966     0     0     0
    72  10060 26081 33450  7590  3612     0     0     0
    73  73696  7097     0     0     0     0     0     0
    74  80793     0     0     0     0     0     0     0
    75  77528  3008   257     0     0     0     0     0
    76  34905 34766 10606   516     0     0     0     0
    77  54799 21256  4116   622     0     0     0     0
    78  57557 21260  1976     0     0     0     0     0
    79  46966 28430  5397     0     0     0     0     0
    80  46394 30478  3921     0     0     0     0     0
    81  80793     0     0     0     0     0     0     0
    82  72153  8640     0     0     0     0     0     0
    83  57732 21901  1160     0     0     0     0     0
    84  73061  7732     0     0     0     0     0     0
    85  78108  2685     0     0     0     0     0     0
    86  64841 15952     0     0     0     0     0     0
    87  75894  4488   411     0     0     0     0     0
    88  77088  3705     0     0     0     0     0     0
    89  48328 32465     0     0     0     0     0     0
    90  80793     0     0     0     0     0     0     0
    91  36647 30720 11100  1868   458     0     0     0
    92  54513 24626  1654     0     0     0     0     0
    93  73946  6847     0     0     0     0     0     0
    94  69274 10406  1113     0     0     0     0     0
    95  43213 20293 14643  2644     0     0     0     0
    96  65015 15778     0     0     0     0     0     0
    97  75428  5365     0     0     0     0     0     0
    98  65968 13004  1821     0     0     0     0     0
    99  78629  2164     0     0     0     0     0     0
    100 52405 15982  9604  2802     0     0     0     0

Beta Diversity

Principal Components Analysis

(For the desktop, the amount of data we generated is too much, so we need to reduce the amount of data we use. To do so we average over time blocks, here of length 100.)

AveragedAbundance <- result$Abundance %>% data.frame(
) %>% dplyr::mutate(
  time = floor(time/100)*100
) %>% dplyr::group_by(
  time
) %>% dplyr::summarise(
  dplyr::across(.fns = ~ mean(.x))
)

We can perform a PCA and see if are data can be summarised by a small number of dimensions. As we shall see, constraining to the first 25 principal components does not harm the system much.

PCA <- prcomp(AveragedAbundance %>% dplyr::select_if(~ any(. > 0)), 
              center = TRUE, scale. = TRUE, rank. = 25)

There is not actually a lot of dependence within the system it appears. Consider the amount of variance explained by the first six principal components (ordered, as is tradition, by amount of variation explained).

head(summary(PCA)$importance[2, ])

    PC1     PC2     PC3     PC4     PC5     PC6 
0.03749 0.03273 0.03095 0.02751 0.02501 0.02482 

The amount of explained variation is as follows.

sum(summary(PCA)$importance[2, ])

[1] 0.99994

The traditional biplot follows, but despite the seeming presence of patterns, so little of the variance is explained that is probably not worth further examination.

ggplot2::autoplot(PCA, loadings = TRUE, 
                  data = AveragedAbundance %>% dplyr::select_if(~ any(. > 0)),
                  colour = "time")

We can try again with presence-absence data instead.

AveragedPA <- result$Abundance %>% data.frame(
) %>% dplyr::mutate(
  time = floor(time/100)*100
) %>% dplyr::mutate(
  dplyr::across(.cols = !time, .fns = ~ .x > 0)
) %>% dplyr::group_by(
  time
) %>% dplyr::summarise(
  dplyr::across(.fns = ~ mean(.x))
)

PCAPA <- prcomp(AveragedPA %>% dplyr::select_if(~ any(. > 0)), 
              center = TRUE, scale. = TRUE, rank. = 25)

head(summary(PCAPA)$importance[2, ])

    PC1     PC2     PC3     PC4     PC5     PC6 
0.05537 0.04584 0.04483 0.03901 0.03578 0.03410 

So it is a bit better, but not by much.

ggplot2::autoplot(PCAPA, loadings = TRUE, 
                  data = AveragedPA %>% dplyr::select_if(~ any(. > 0)),
                  colour = "time")

I should note that this is not a failure of the method persay. What this says to me is that dimension reduction of the system cannot reduce the system to a human-readable set of descriptors. The system is still being reduced (from 100 Species * 10 Environments to 25 Components to cover

sum(summary(PCA)$importance[2, ])

[1] 0.99994

of the variation). My initial hypothesis for when the systems are coupled is that, as coupling increases, the number of principal components should decrease, since one system’s changes will be better predictors of the next system’s changes. (An interesting related question; do the number of principal components correlate with the amount of biodiversity in the system?

ifrm(AveragedAbundance)
ifrm(AveragedPA)
ifrm(PCA)
ifrm(PCAPA)

Beta Diversity Over Time

Since trying to reduce the system dimensionally does not work (well enough) a next attempt might be to consider how the pair-wise beta diversity changes over the course of the simulation. Initial problems include that there are quite a few measures of beta diversity as well as that the (square of the) number of environments will determine how many entries we need to consider.

To try this, we are going to reorganise our data into data frame with attributes of environment, time, and traits/species.

Environments <- lapply(
  1:result$NumEnvironments,
  function(i, abund, numSpecies) {
    time <- abund[, 1]
    env <- abund[, 1 + 1:numSpecies + numSpecies * (i - 1)]
    
    retval <- data.frame(Environment = toString(i),
                         Time = time
    )
    
    retval <- cbind(retval, env)
    colnames(retval) <- c("Environment", "Time", 
                          paste0("Basal", 1:34), 
                          paste0("Consumer", 35:100)
    )
    return(retval)
                          
  },
  abund = result$Abundance %>% data.frame(
  ) %>% dplyr::mutate(
    time = floor(time/100)*100
  ) %>% dplyr::group_by(
    time
  ) %>% dplyr::summarise(
    dplyr::across(.fns = ~ mean(.x))
  ),
  numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)

Environments <- dplyr::bind_rows(Environments)
Environments <- Environments %>% dplyr::filter(
  dplyr::if_any(dplyr::starts_with(c("Basal", "Consumer")), ~ (.x != 0)),
)

Cluster Analysis Over Time

One perhaps interesting idea is to try to group the environments by cluster by considering the abundances (or presence-absence) of the species present as traits. The question then is whether the environments have a tendency to attract to specific points or if they wander around each other without relation. (The latter is more neutral.) If they appear to be convergent (which so far would seem to disagree mostly with the alpha diversity analyses), then that would imply that dynamics determine the majority of the system’s fate, while if they instead wander more randomly, then they would appear to be dominated by the neutral mechanisms. (Of course, this is affected by the rate of the neutral mechanisms, so it exists on a definite continuum.) (Note, of course, that cluster analysis usually includes something like PCA, so we would not necessarily expect a different result.)

My working hypothesis for how to use the distance measures is that we need a metric measure (Kindt and Coe, 2005) with any identification scrubbed off (site, time collected). Kindt and Coe suggest that taking the square root of Bray-Curtis makes it metric and thus usable for mathematical distance analyses (such as hierarchical clustering). The vegan package recomends using Jaccard instead of Bray-Curtis, for the same reason.

EnvironmentDistance <- Environments %>% dplyr::select(
  -Environment, -Time
) %>% vegan::vegdist(method = "jaccard")

EnvDistClust <- hclust(EnvironmentDistance)

plot(EnvDistClust, labels = FALSE)

Note that are major breakdowns extremely early yielding a large number of clusters but there are expected to be 10 environments (if historical contingency mattered most) or a small number of clusters (if they are converging). The large number of clusters says to me that the neutral dynamics are jumping rapidly from cluster to cluster.

We can try again with presence absence instead.

EnvironmentDistancePA <- Environments %>% dplyr::select(
  -Environment, -Time
) %>% vegan::vegdist(method = "jaccard", binary = TRUE)

EnvDistClustPA <- hclust(EnvironmentDistancePA)

plot(EnvDistClustPA, labels = FALSE)

Results, Line Spatial Structure

ifrm(Diversity)
ifrm(EnvDistClust)
ifrm(EnvDistClustPA)
ifrm(EnvDiversity)
ifrm(Environments)
ifrm(EnvironmentDistance)
ifrm(EnvironmentDistancePA)

We repeat the procedures above for the linear set. We use the same pool and interaction matrices, but change the space so that ``adjacent’’ communities can disperse (e.g. 1 <-> 2 <-> 3 <->…<-> 9 <-> 10). Similarly, the history is the same in terms of arrivals and extinctions.

load(file.path(
  "..", "experiments", "MNA-FirstAttempt-Result-Env10-Line.RData")
)
toRemove <- result$Abundance <= result$Parameters$EliminationThreshold
result$Abundance[toRemove] <- 0
ifrm(toRemove)

Abundance

With 10 environments, it is probably not helpful to check 10 individual abundance curves, but looking at the first one might be helpful.

LawMorton1996_PlotAbundance(result$Abundance[seq(from = 1, 
                                                 to = nrow(result$Abundance), 
                                                 by = 10), c(1, 2:101)]) -> obj;

obj + ggplot2::scale_y_log10() + ggplot2::guides(color = FALSE)

Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
Warning: Transformation introduced infinite values in continuous y-axis

Every vertical line is a species introduction or extinction by neutral dynamics. Note that, due to dispersal, these should lose some of their abundance to spreading, but if that spreading is faster than the population can recoup its losses, than it can invade and be wiped out (since invadability does not consider spatial dynamics, only local ecological dynamics). The primary notes then are that there are a large number more events appearing to occur (due to the spatial arrangement), but the system as a whole is stabler (since the species that get knocked out can be recaptured from space as well).

Alpha Diversity plots

Intro

Diversity <- lapply(
  1:result$NumEnvironments,
  function(i, abund, numSpecies) {
    time <- abund[, 1]
    env <- abund[, 1 + 1:numSpecies + numSpecies * (i - 1)]
    richness <- rowSums(env != 0)
    abundSum <- rowSums(env)
    entropy <- env / abundSum
    entropy <- - apply(
      entropy, MARGIN = 1,
      FUN = function(x) {
        sum(ifelse(x != 0, x * log(x), 0))
      })
    species <- apply(
      env, MARGIN = 1,
      FUN = function(x) {
        toString(which(x > 0))
      }
    )
    data.frame(Time = time, 
               Richness = richness, 
               Entropy = entropy,
               Species = species,
               Environment = i)
  },
  abund = result$Abundance,
  numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)


Diversity <- dplyr::bind_rows(Diversity)
Diversity <- Diversity %>% dplyr::mutate(
  Evenness = Entropy / log(Richness)
)

Richness

ggplot2::ggplot(
  Diversity, 
  ggplot2::aes(
    x = Time,
    y = Richness,
    color = factor(Environment),
    alpha = ifelse(Environment == 1, 1, 0.3)
  )
) + ggplot2::geom_line(
) + ggplot2::geom_line(
  data = Diversity %>% dplyr::group_by(
    Time
    ) %>% dplyr::summarise(
      Richness = mean(Richness)
    ),
  mapping = ggplot2::aes(
    x = Time,
    y = Richness
  ),
  color = "black",
  inherit.aes = FALSE
) + ggplot2::guides(
  alpha = "none"
) + ggplot2::scale_color_discrete(
  "Environment"
)

Entropy

ggplot2::ggplot(
  Diversity, 
  ggplot2::aes(
    x = Time,
    y = Entropy,
    color = factor(Environment),
    alpha = ifelse(Environment == 1, 1, 0.3)
  )
) + ggplot2::geom_line(
) + ggplot2::geom_line(
  data = Diversity %>% dplyr::group_by(
    Time
    ) %>% dplyr::summarise(
      Entropy = mean(Entropy)
    ),
  mapping = ggplot2::aes(
    x = Time,
    y = Entropy
  ),
  color = "black",
  inherit.aes = FALSE
) + ggplot2::guides(
  alpha = "none"
) + ggplot2::scale_color_discrete(
  "Environment"
)

Warning: Removed 1123 row(s) containing missing values (geom_path).
Warning: Removed 110 row(s) containing missing values (geom_path).

Entropy-Richness

plotly::plot_ly(data = Diversity %>% dplyr::filter(Environment < 2),
                x = ~Richness, y = ~Entropy, z = ~Time, type = "scatter3d",
                mode = "lines", opacity = 1, line = list(color = ~Time))

Evenness

ggplot2::ggplot(
  Diversity, 
  ggplot2::aes(
    x = Time,
    y = Evenness,
    color = factor(Environment),
    alpha = ifelse(Environment == 1, 1, 0.3)
  )
) + ggplot2::geom_line(
) + ggplot2::geom_line(
  data = Diversity %>% dplyr::group_by(
    Time
    ) %>% dplyr::summarise(
      Evenness = mean(Evenness)
    ),
  mapping = ggplot2::aes(
    x = Time,
    y = Evenness
  ),
  color = "black",
  inherit.aes = FALSE
) + ggplot2::guides(
  alpha = "none"
) + ggplot2::scale_color_discrete(
  "Environment"
)

Warning: Removed 1440 row(s) containing missing values (geom_path).
Warning: Removed 136 row(s) containing missing values (geom_path).

Environment Diversity

EnvDiversity <- lapply(
    1:((ncol(result$Abundance) - 1) / result$NumEnvironments),
    function(i, abund, numSpecies) {
        time <- abund[, 1]
        env <- abund[, 1 + i + numSpecies * (1:result$NumEnvironments - 1)]
        richness <- rowSums(env != 0)
        abundSum <- rowSums(env)
        environments <- apply(
            env, MARGIN = 1,
            FUN = function(x) {
                toString(which(x > 0))
            }
        )
        data.frame(Time = time, 
                   Richness = richness, 
                   Abundance = abundSum,
                   Species = i,
                   Environments = environments)
    },
    abund = result$Abundance,
    numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)

EnvDiversity <- dplyr::bind_rows(EnvDiversity)

ggplot2::ggplot(
  EnvDiversity %>% dplyr::filter(Richness > 1), 
  aes(x = Time, y = Richness, color = Species)
) + geom_point(
  alpha = 0.01, size = 3
) + guides(
  color = "none"
)

with(EnvDiversity %>% mutate(
  Time = floor(Time)
  ) %>% group_by(
    Time, Species
    ) %>% summarise(
      Richness = round(mean(Richness)),
      .groups = "drop"
      ),
     table(Species, Richness))

       Richness
Species     0     1     2     3     4     5     6     7     8     9    10
    1   80793     0     0     0     0     0     0     0     0     0     0
    2   80624     0     1     1     2     1     1     2     0     2   159
    3   80793     0     0     0     0     0     0     0     0     0     0
    4   80047     0     0     2     0     1     0     0     1     0   742
    5   79905     0     6     0     6     2     5     1     0    11   857
    6   80793     0     0     0     0     0     0     0     0     0     0
    7   80793     0     0     0     0     0     0     0     0     0     0
    8   80793     0     0     0     0     0     0     0     0     0     0
    9   80793     0     0     0     0     0     0     0     0     0     0
    10  80793     0     0     0     0     0     0     0     0     0     0
    11  80793     0     0     0     0     0     0     0     0     0     0
    12    680     0     0     0     0     0     0     0     0     1 80112
    13  80639     1     1     1     1     1     1     2     1     1   144
    14     84     0     0     0     0     0     0     0     0     0 80709
    15  80793     0     0     0     0     0     0     0     0     0     0
    16  80793     0     0     0     0     0     0     0     0     0     0
    17  80793     0     0     0     0     0     0     0     0     0     0
    18  80793     0     0     0     0     0     0     0     0     0     0
    19  80793     0     0     0     0     0     0     0     0     0     0
    20  80793     0     0     0     0     0     0     0     0     0     0
    21  80793     0     0     0     0     0     0     0     0     0     0
    22     62     0     0     0     0     0     0     0     0     1 80730
    23  80793     0     0     0     0     0     0     0     0     0     0
    24  80546     0     1     0     0     0     0     1     0     0   245
    25  80793     0     0     0     0     0     0     0     0     0     0
    26  80793     0     0     0     0     0     0     0     0     0     0
    27  80793     0     0     0     0     0     0     0     0     0     0
    28  80793     0     0     0     0     0     0     0     0     0     0
    29  80721     1     0     0     1     0     0     0     0     1    69
    30  69391     9     2     8     1     8     5     2     8    12 11347
    31  80793     0     0     0     0     0     0     0     0     0     0
    32    500     0     0     0     0     0     1     0     0     1 80291
    33  66491     5    15    10    14     2    14     5    14    17 14206
    34  80793     0     0     0     0     0     0     0     0     0     0
    35  80332     1     0     1     1     1     1     1     0     2   453
    36  80793     0     0     0     0     0     0     0     0     0     0
    37  46279     1     1     2     1     1     3     1     2     4 34498
    38  80793     0     0     0     0     0     0     0     0     0     0
    39  80793     0     0     0     0     0     0     0     0     0     0
    40  80793     0     0     0     0     0     0     0     0     0     0
    41  80793     0     0     0     0     0     0     0     0     0     0
    42  80793     0     0     0     0     0     0     0     0     0     0
    43  27630     2     2     0     4     0     0     1     4     7 53143
    44  60232     1     8     3     3     7     2     7     5    14 20511
    45  80793     0     0     0     0     0     0     0     0     0     0
    46  80793     0     0     0     0     0     0     0     0     0     0
    47  52191     5     8     3     7     3     7     5     9    12 28543
    48  80793     0     0     0     0     0     0     0     0     0     0
    49  80793     0     0     0     0     0     0     0     0     0     0
    50  80793     0     0     0     0     0     0     0     0     0     0
    51  80793     0     0     0     0     0     0     0     0     0     0
    52  80793     0     0     0     0     0     0     0     0     0     0
    53  80793     0     0     0     0     0     0     0     0     0     0
    54  80793     0     0     0     0     0     0     0     0     0     0
    55  78735     1     1     2     0     1     2     0     0     2  2049
    56  80793     0     0     0     0     0     0     0     0     0     0
    57  79020     0     0     1     0     0     0     0     1     0  1771
    58  80793     0     0     0     0     0     0     0     0     0     0
    59  80793     0     0     0     0     0     0     0     0     0     0
    60  80793     0     0     0     0     0     0     0     0     0     0
    61  80466     0     1     0     0     1     1     1     0     1   322
    62  80793     0     0     0     0     0     0     0     0     0     0
    63  80793     0     0     0     0     0     0     0     0     0     0
    64  80793     0     0     0     0     0     0     0     0     0     0
    65  80793     0     0     0     0     0     0     0     0     0     0
    66  80793     0     0     0     0     0     0     0     0     0     0
    67  80793     0     0     0     0     0     0     0     0     0     0
    68  79023     1     5     4     7     5     6     4     4     4  1730
    69  78135     2     2     1     1     1     6     1     2     5  2637
    70  80793     0     0     0     0     0     0     0     0     0     0
    71   1394     0     0     0     0     0     0     0     0     1 79398
    72    292     1     0     0     0     0     1     0     1     0 80498
    73  80793     0     0     0     0     0     0     0     0     0     0
    74  80793     0     0     0     0     0     0     0     0     0     0
    75  80793     0     0     0     0     0     0     0     0     0     0
    76  78676     1     4     1     4     0     6     1     3     6  2091
    77  79765     0     0     1     0     0     1     0     0     1  1025
    78  79253     0     0     1     0     1     0     0     0     0  1538
    79  80793     0     0     0     0     0     0     0     0     0     0
    80  80793     0     0     0     0     0     0     0     0     0     0
    81  80793     0     0     0     0     0     0     0     0     0     0
    82  80793     0     0     0     0     0     0     0     0     0     0
    83  80793     0     0     0     0     0     0     0     0     0     0
    84  80793     0     0     0     0     0     0     0     0     0     0
    85  80793     0     0     0     0     0     0     0     0     0     0
    86  64862     3     6     4     7     4     6     3     8     6 15884
    87  80793     0     0     0     0     0     0     0     0     0     0
    88  80793     0     0     0     0     0     0     0     0     0     0
    89  80793     0     0     0     0     0     0     0     0     0     0
    90  80793     0     0     0     0     0     0     0     0     0     0
 [ reached getOption("max.print") -- omitted 10 rows ]

Beta Diversity

Principal Components Analysis

AveragedAbundance <- result$Abundance %>% data.frame(
) %>% dplyr::mutate(
  time = floor(time/100)*100
) %>% dplyr::group_by(
  time
) %>% dplyr::summarise(
  dplyr::across(.fns = ~ mean(.x))
)

PCA <- prcomp(AveragedAbundance %>% dplyr::select_if(~ any(. > 0)), 
              center = TRUE, scale. = TRUE, rank. = 25)

head(summary(PCA)$importance[2, ])

    PC1     PC2     PC3     PC4     PC5     PC6 
0.19145 0.11089 0.10844 0.07943 0.05993 0.04530 

sum(summary(PCA)$importance[2, ])

[1] 0.99996

ggplot2::autoplot(PCA, loadings = TRUE, 
                  data = AveragedAbundance %>% dplyr::select_if(~ any(. > 0)),
                  colour = "time")

AveragedPA <- result$Abundance %>% data.frame(
) %>% dplyr::mutate(
  time = floor(time/100)*100
) %>% dplyr::mutate(
  dplyr::across(.cols = !time, .fns = ~ .x > 0)
) %>% dplyr::group_by(
  time
) %>% dplyr::summarise(
  dplyr::across(.fns = ~ mean(.x))
)

PCAPA <- prcomp(AveragedPA %>% dplyr::select_if(~ any(. > 0)), 
              center = TRUE, scale. = TRUE, rank. = 25)

head(summary(PCAPA)$importance[2, ])

    PC1     PC2     PC3     PC4     PC5     PC6 
0.21750 0.12185 0.08092 0.05943 0.05215 0.04468 

ggplot2::autoplot(PCAPA, loadings = TRUE, 
                  data = AveragedPA %>% dplyr::select_if(~ any(. > 0)),
                  colour = "time")

sum(summary(PCA)$importance[2, ])

[1] 0.99996

ifrm(AveragedAbundance)
ifrm(AveragedPA)
ifrm(PCA)
ifrm(PCAPA)

Beta Diversity Over Time

Environments <- lapply(
  1:result$NumEnvironments,
  function(i, abund, numSpecies) {
    time <- abund[, 1]
    env <- abund[, 1 + 1:numSpecies + numSpecies * (i - 1)]
    
    retval <- data.frame(Environment = toString(i),
                         Time = time
    )
    
    retval <- cbind(retval, env)
    colnames(retval) <- c("Environment", "Time", 
                          paste0("Basal", 1:34), 
                          paste0("Consumer", 35:100)
    )
    return(retval)
                          
  },
  abund = result$Abundance %>% data.frame(
  ) %>% dplyr::mutate(
    time = floor(time/100)*100
  ) %>% dplyr::group_by(
    time
  ) %>% dplyr::summarise(
    dplyr::across(.fns = ~ mean(.x))
  ),
  numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)

Environments <- dplyr::bind_rows(Environments)
Environments <- Environments %>% dplyr::filter(
  dplyr::if_any(dplyr::starts_with(c("Basal", "Consumer")), ~ (.x != 0)),
)

Cluster Analysis Over Time

EnvironmentDistance <- Environments %>% dplyr::select(
  -Environment, -Time
) %>% vegan::vegdist(method = "jaccard")

EnvDistClust <- hclust(EnvironmentDistance)

plot(EnvDistClust, labels = FALSE)

EnvironmentDistancePA <- Environments %>% dplyr::select(
  -Environment, -Time
) %>% vegan::vegdist(method = "jaccard", binary = TRUE)

EnvDistClustPA <- hclust(EnvironmentDistancePA)

plot(EnvDistClustPA, labels = FALSE)

Results, Slow Line Spatial Structure

ifrm(Diversity)
ifrm(EnvDistClust)
ifrm(EnvDistClustPA)
ifrm(EnvDiversity)
ifrm(Environments)
ifrm(EnvironmentDistance)
ifrm(EnvironmentDistancePA)

We repeat the procedures above for the linear set. We use the same pool and interaction matrices, but change the space so that ``adjacent’’ communities can disperse (e.g. 1 <-> 2 <-> 3 <->…<-> 9 <-> 10). Similarly, the history is the same in terms of arrivals and extinctions.

load(file.path(
  "..", "experiments", "MNA-FirstAttempt1e+06-Result-Env10-Line.RData")
)
toRemove <- result$Abundance <= result$Parameters$EliminationThreshold
result$Abundance[toRemove[1:(length(toRemove)/2)]] <- 0
result$Abundance[toRemove[(length(toRemove)/2 + 1):length(toRemove)]] <- 0
ifrm(toRemove)

Abundance

With 10 environments, it is probably not helpful to check 10 individual abundance curves, but looking at the first one might be helpful.

LawMorton1996_PlotAbundance(result$Abundance[seq(from = 1, 
                                                 to = nrow(result$Abundance), 
                                                 by = 10), c(1, 2:101)]) -> obj;

obj + ggplot2::scale_y_log10() + ggplot2::guides(color = FALSE)

Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
Warning: Transformation introduced infinite values in continuous y-axis

Alpha Diversity plots

Intro

Diversity <- lapply(
  1:result$NumEnvironments,
  function(i, abund, numSpecies) {
    time <- abund[, 1]
    env <- abund[, 1 + 1:numSpecies + numSpecies * (i - 1)]
    richness <- rowSums(env != 0)
    abundSum <- rowSums(env)
    entropy <- env / abundSum
    entropy <- - apply(
      entropy, MARGIN = 1,
      FUN = function(x) {
        sum(ifelse(x != 0, x * log(x), 0))
      })
    species <- apply(
      env, MARGIN = 1,
      FUN = function(x) {
        toString(which(x > 0))
      }
    )
    data.frame(Time = time, 
               Richness = richness, 
               Entropy = entropy,
               Species = species,
               Environment = i)
  },
  abund = result$Abundance,
  numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)


Diversity <- dplyr::bind_rows(Diversity)
Diversity <- Diversity %>% dplyr::mutate(
  Evenness = Entropy / log(Richness)
)

Richness

ggplot2::ggplot(
  Diversity, 
  ggplot2::aes(
    x = Time,
    y = Richness,
    color = factor(Environment),
    alpha = ifelse(Environment == 1, 1, 0.3)
  )
) + ggplot2::geom_line(
) + ggplot2::geom_line(
  data = Diversity %>% dplyr::group_by(
    Time
    ) %>% dplyr::summarise(
      Richness = mean(Richness)
    ),
  mapping = ggplot2::aes(
    x = Time,
    y = Richness
  ),
  color = "black",
  inherit.aes = FALSE
) + ggplot2::guides(
  alpha = "none"
) + ggplot2::scale_color_discrete(
  "Environment"
)

geom_path: Each group consists of only one observation. Do you need to adjust the group
aesthetic?

Entropy

ggplot2::ggplot(
  Diversity, 
  ggplot2::aes(
    x = Time,
    y = Entropy,
    color = factor(Environment),
    alpha = ifelse(Environment == 1, 1, 0.3)
  )
) + ggplot2::geom_line(
) + ggplot2::geom_line(
  data = Diversity %>% dplyr::group_by(
    Time
    ) %>% dplyr::summarise(
      Entropy = mean(Entropy)
    ),
  mapping = ggplot2::aes(
    x = Time,
    y = Entropy
  ),
  color = "black",
  inherit.aes = FALSE
) + ggplot2::guides(
  alpha = "none"
) + ggplot2::scale_color_discrete(
  "Environment"
)

Warning: Removed 197949 row(s) containing missing values (geom_path).
Warning: Removed 1 row(s) containing missing values (geom_path).

Entropy-Richness

plotly::plot_ly(data = Diversity %>% dplyr::filter(Environment < 2),
                x = ~Richness, y = ~Entropy, z = ~Time, type = "scatter3d",
                mode = "lines", opacity = 1, line = list(color = ~Time))

Evenness

ggplot2::ggplot(
  Diversity, 
  ggplot2::aes(
    x = Time,
    y = Evenness,
    color = factor(Environment),
    alpha = ifelse(Environment == 1, 1, 0.3)
  )
) + ggplot2::geom_line(
) + ggplot2::geom_line(
  data = Diversity %>% dplyr::group_by(
    Time
    ) %>% dplyr::summarise(
      Evenness = mean(Evenness)
    ),
  mapping = ggplot2::aes(
    x = Time,
    y = Evenness
  ),
  color = "black",
  inherit.aes = FALSE
) + ggplot2::guides(
  alpha = "none"
) + ggplot2::scale_color_discrete(
  "Environment"
)

Warning: Removed 294140 row(s) containing missing values (geom_path).
Warning: Removed 1 row(s) containing missing values (geom_path).

Environment Diversity

EnvDiversity <- lapply(
    1:((ncol(result$Abundance) - 1) / result$NumEnvironments),
    function(i, abund, numSpecies) {
        time <- abund[, 1]
        env <- abund[, 1 + i + numSpecies * (1:result$NumEnvironments - 1)]
        richness <- rowSums(env != 0)
        abundSum <- rowSums(env)
        environments <- apply(
            env, MARGIN = 1,
            FUN = function(x) {
                toString(which(x > 0))
            }
        )
        data.frame(Time = time, 
                   Richness = richness, 
                   Abundance = abundSum,
                   Species = i,
                   Environments = environments)
    },
    abund = result$Abundance,
    numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)

EnvDiversity <- dplyr::bind_rows(EnvDiversity)

ggplot2::ggplot(
  EnvDiversity %>% dplyr::filter(Richness > 1), 
  aes(x = Time, y = Richness, color = Species)
) + geom_point(
  alpha = 0.01, size = 3
) + guides(
  color = "none"
)

with(EnvDiversity %>% mutate(
  Time = floor(Time)
  ) %>% group_by(
    Time, Species
    ) %>% summarise(
      Richness = round(mean(Richness)),
      .groups = "drop"
      ),
     table(Species, Richness))

       Richness
Species 0 1 4 5
    1   1 0 0 0
    2   1 0 0 0
    3   1 0 0 0
    4   1 0 0 0
    5   1 0 0 0
    6   1 0 0 0
    7   1 0 0 0
    8   1 0 0 0
    9   1 0 0 0
    10  1 0 0 0
    11  1 0 0 0
    12  0 0 0 1
    13  1 0 0 0
    14  0 0 1 0
    15  1 0 0 0
    16  1 0 0 0
    17  1 0 0 0
    18  1 0 0 0
    19  1 0 0 0
    20  1 0 0 0
    21  1 0 0 0
    22  0 0 0 1
    23  1 0 0 0
    24  1 0 0 0
    25  1 0 0 0
    26  1 0 0 0
    27  1 0 0 0
    28  1 0 0 0
    29  1 0 0 0
    30  1 0 0 0
    31  1 0 0 0
    32  0 1 0 0
    33  1 0 0 0
    34  1 0 0 0
    35  1 0 0 0
    36  1 0 0 0
    37  1 0 0 0
    38  1 0 0 0
    39  1 0 0 0
    40  1 0 0 0
    41  1 0 0 0
    42  1 0 0 0
    43  1 0 0 0
    44  1 0 0 0
    45  1 0 0 0
    46  1 0 0 0
    47  0 1 0 0
    48  1 0 0 0
    49  1 0 0 0
    50  1 0 0 0
    51  1 0 0 0
    52  1 0 0 0
    53  1 0 0 0
    54  1 0 0 0
    55  1 0 0 0
    56  1 0 0 0
    57  1 0 0 0
    58  1 0 0 0
    59  1 0 0 0
    60  1 0 0 0
    61  1 0 0 0
    62  1 0 0 0
    63  1 0 0 0
    64  1 0 0 0
    65  1 0 0 0
    66  1 0 0 0
    67  1 0 0 0
    68  1 0 0 0
    69  1 0 0 0
    70  1 0 0 0
    71  1 0 0 0
    72  0 1 0 0
    73  1 0 0 0
    74  1 0 0 0
    75  1 0 0 0
    76  1 0 0 0
    77  1 0 0 0
    78  1 0 0 0
    79  1 0 0 0
    80  1 0 0 0
    81  1 0 0 0
    82  1 0 0 0
    83  1 0 0 0
    84  1 0 0 0
    85  1 0 0 0
    86  1 0 0 0
    87  1 0 0 0
    88  1 0 0 0
    89  1 0 0 0
    90  1 0 0 0
    91  1 0 0 0
    92  1 0 0 0
    93  1 0 0 0
    94  1 0 0 0
    95  1 0 0 0
    96  1 0 0 0
    97  1 0 0 0
    98  1 0 0 0
    99  1 0 0 0
    100 1 0 0 0

Beta Diversity

Principal Components Analysis

AveragedAbundance <- result$Abundance %>% data.frame(
) %>% dplyr::mutate(
  time = floor(time/100)*100
) %>% dplyr::group_by(
  time
) %>% dplyr::summarise(
  dplyr::across(.fns = ~ mean(.x))
)

PCA <- prcomp(AveragedAbundance %>% dplyr::select_if(~ any(. > 0)), 
              center = TRUE, scale. = TRUE, rank. = 25)

Error in prcomp.default(AveragedAbundance %>% dplyr::select_if(~any(. >  : 
  cannot rescale a constant/zero column to unit variance

Beta Diversity Over Time

Cluster Analysis Over Time